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OPTIMAL PERTURBATIONS OF RUNGE-KUTTA METHODS 

INMACULADA HIGUERAS*, DAVID I. KETCHESONt, AND TIHAMER A. KOCSIS* 


Abstract. Perturbed Runge-Kutta methods (also referred to as downwind Runge—Kutta methods) can guaran- 
tee monotonicity preservation under larger step sizes relative to their traditional Runge-Kutta counterparts. In this 
paper we study, the question of how to optimally perturb a given method in order to increase the radius of absolute 
monotonicity (a.m.). We prove that for methods with zero radius of a.m., it is always possible to give a perturbation 
with positive radius. We first study methods for linear problems and then methods for nonlinear problems. In each 
case, we prove upper bounds on the radius of a.m., and provide algorithms to compute optimal perturbations. We 
also provide optimal perturbations for many known methods. 
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1. Introduction. Strong stability preserving Runge-Kutta (RK) methods were first intro¬ 
duced by Shu and Osher |23| in the context of time integration for first-order hyperbolic conservation 
laws: 


Ut + F{U)^ = Q, U{x,t = Q)=Uo. (1.1) 

After semi-discretization, 0 takes the form of an initial-value ordinary differential equation 
system: 


u'{t) = f{u) u(0) = uo , (1.2) 

where / is a discrete approximation to —J-. In the scalar case lA is dissipative, and it is natural to 
seek a semi-discretization that is dissipative: 


|hll<0, (1.3) 

where || • || denotes a convex functional (e.g., a norm, a semi-norm, ...). This is achieved by biasing 
the discretization / in the upwind direction. A necessary condition for ( |1.3[ ) is monotonicity under 
an explicit Euler step [T^ p. 501]: 

Iju -I- hf{v)\\ < ||f;||, for all u, and for h satisfying 0 < h < ho, (1.4) 

where ho > 0 (in general ho may depend on v). Let denote approximations, computed by 

some numerical integrator, to the solution at successive time steps and tn+i = tn + h. Under 
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the forward Euler monotonicity condition (1.4), it is possible to prove that many Runge-Kutta and 
linear multistep methods also give monotone solutions; i.e., solutions that satisfy 


||m„+i|| < ||m„ 


for h satisfying 0 < h < Rho. 


(1.5) 


Such methods are known as strong stability preserving (SSP) methods, and the factor R is known 
as the radius of absolute monotonicity or SSP coefficient of the method. SSP methods necessarily 


have non-negative coefficients, since the monotonicity property is proved using (1.4) and convexity. 


Monotonicity cannot be ensured using only assumption (1.4) in general for methods with neg¬ 


ative coefficients Thm. 4.2], or even for some methods (such as the classical fourth-order RK 
method) with non-negative coefficients [T^ Thm. 9.6]. In order to accommodate such methods, a 
second discrete approximation to —J- is introduced and referred to as /. This discretization must 
be monotone under an explicit Euler step with negative step size: 


\\v-hf{v)\\ < Hull 


for all V, and for h satisfying 0 < h < ho, 


( 1 . 6 ) 


where ho > 0. In the context of hyperbolic problems, / must be biased in the downwind direction, 
and typically = h^. The downwind spatial discretization / is to be used in place of / wherever a 
negative coefficient appears in the time integration method, in order to ensure monotonicity of the 
overall method. Introduction of a downwind discretization makes it possible to ensure monotonicity 
for a broader class of methods, including the classical RK method of order four. It also makes it 
possible to ensure monotonicity for many methods under larger step sizes. 

Methods that use both upwind and downwind operators can naturally be viewed as perturbed 
Runge-Kutta methods. Although they are also connected to additive RK methods (see Einni), 
in the present work we will employ the perturbation viewpoint, and refer to methods that use 
downwind discretization as perturbed RK methods. 

During the last quarter century, a number of additional authors have studied monotonicity 
for methods that use downwind discretization. The main motivation for this work has been to 
break the “order barrier” that restricts explicit RK methods to order four, or to find new methods 
with larger SSP coefficients. In this context, numerical optimization of the SSP coefficient for RK 
methods with negative coefficients was conducted for explicit methods in mmiis] and for implicit 
methods in m- In each case, optimization was carried out over methods with a specified order 
and number of stages. 

The present work stems from a different motivation. Monotonicity preservation is not the 
only numerical property of interest in applications, and practitioners may wish to use a particular 
integrator that has small or zero SSP coefficient. Our goal is then to perturb the prescribed 
method in order to achieve larger monotonicity-preserving timesteps. Little work has been done in 
this direction, because it is not known how to find the best perturbation for a given method. That 
problem is the main focus of this work. Most of our results concern only explicit methods, although 
one major result (Theorem |3.1[) pertains also to implicit methods. 


1.1. Perturbed Runge-Kutta methods. A Runge-Kutta method applied to the initial 
value problem (1.2) computes approximations Un « u{tn) by 


(1.7a) 
(1.7b) 
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Y = Une -f hKF , 

Un+l = Kjj+i. 






Here s is the number of stages, e is a vector whose entries are equal to one, Y is the vector containing 
the stage values and the numerical solution, Y = (Fi, ..., Fg+i)*, [F]i = /(F^), and K is the 

(s + 1) X (s + 1) matrix of Butcher coefficients: 


K = 




In this work we study perturbed Runge-Kutta methods: 

(1.8a) 
(1.8b) 


F = M„e + hKF + hK{F - F ), 

^n+l ~ ^+1; 


where K is given by 



and F is deffned analogously to F. We assume that the perturbation A has the same structure 
(strictly lower-triangular, lower-triangular, or full) as the matrix A. 

Observe that the perturbed method (1.8) reduces to the RK method (1.71 when F = F. Method 


(1.8) may be viewed as approximating the solution of the perturbed problem 

u'it) = f{u) + {f{u) - f{u)), 


where f ~ f ■ 

We assume that / and / satisfy the explicit Euler assumptions (1.4) and (1.6), respectively, 
with ho = Hq. 


Most previous works, including ElIS], have focused on methods with the following property 
Property C: ^ 

We say that a perturbation K to an RK method K possesses property C if, for 
each value of j 


Kij A 0 (for some i) 


Kij = 0 (for all i). 


(1.9) 


In words, property C means that in the _)th column, only one of K, K has any nonzero entries. Thus, 
only one of f{yj),fiyj) need ever be evaluated, so only s total function evaluations are required 
per step. In [6] it was shown that for WENO discretizations, the cost of computing both f{yj) and 
fipj) is much less than twice the cost of computing f{yj) alone. Therefore methods that without 
property C may also be of practical interest. In the present work, we do not assume property C. 


1.2. Scope and outline. The central question of the present work is 
• Given a ffxed RK method, what perturbation results in the largest radius of absolute 
monotonicity for the perturbed method? 

We investigate this question in the context of both linear problems (Section and nonlinear 
problems (Section]^. 

For explicit methods applied to linear problems, the question above can be cast in terms of 
absolute monotonicity of the (bivariate) stability polynomial. In Section 2.1.1 we prove a general 
upper bound on the radius of absolute monotonicity of the stability polynomial of an explicit 
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perturbed RK method with s stages and linear order p. In Section |2.1.2 we provide an algorithm 
for computing tighter bounds, and tabulate some of the resulting numerical values. Examples of 
optimal methods are given in |2.2[ 

Section]^ is devoted to perturbed Runge-Kutta methods for nonlinear problems. Theorem |3.1| 
states that a perturbation with positive radius of absolute monotonicity exists for every Runge- 
Kutta method. Theorems |3.3| and |3.4| g ive simple upper bounds on the optimal perturbed radius for 
a given explicit method. In Section f3.5| we give two algorithms for computing optimal perturbations. 
The first is provably correct but approximate, while the second is heuristic but exact and agrees 
with the first in all cases we have tested. Both are applicable only to explicit methods. These 
algorithms have been implemented in the software package Nodepy [m. We conclude Section d 
with an application of the theorems and algorithms to optimal perturbations of some RK methods 
from the literature. Among the results is the first truly optimal perturbation for the classical 
4th-order method of Kutta. 

Section [^contains some conclusions as well as some open questions to be studied in the future. 

Finally, in the Appendix we give some details on perturbations for the family of second order 
2-stage methods and the classical fourth order RK method. 

2. Explicit perturbed Runge Kutta methods for linear problems. To study the be¬ 


havior of the perturbed Runge-Kutta method (1.8) for linear problems, we apply it to a linear 
scalar test problem, setting f{u) = Xu and f{u) = Xu in (|1.8[). This results in the iteration 


^-|-l — z)Un, 


where z = hX, z = hX and 

z) = 1 + -I- (z -I- z)6‘^ (^I - zA - {z + z)A^ 


-1 


e. 


( 2 . 1 ) 


We refer to (2.11 as the stability function of the perturbed Runge-Kutta method (1.8). 


We say that a function i/; : M —>■ K is absolutely monotonic (a.m.) at ^ if all derivatives at ^ 
exist and they are non-negative. For a function ■)/) : —>■ K the concept of absolute monotonicity 

can be defined in a similar way: '0 is a.m. at (^,^) if all derivatives at (5,^) exist and they are 
non-negative. 

Given a function 'ip[z,z), we define the radius of absolute monotonicity as 


= sup {r S K I r = 0, or r > 0 , and V’(z, z) is a.m. at (—r, —r)} . 


( 2 . 2 ) 


For a perturbed Runge-Kutta method (1.8) with coefficients (AT, AT), we write K) to 


denote the radius of absolute monotonicity of its stability function: 

RuniK,K) = if)) . 

The quantity i?Lin(-f^)^) is referred to as the threshold factor due to its role in the step size for 
monotonicity. For a given Runge-Kutta method (1.7) with coefficients A', we are interested in 


determining perturbations K that give the largest threshold factor. The corresponding threshold 
factor of the optimal perturbation is denoted by 


^Un(^)=SUpi?Lin(^,^). 

K 


(2.3) 
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The supremum in (2.31 is taken over all strictly lower triangular matrices K, in order to preserve 
the explicit nature of the method. A perturbation K such that 


Rl^^{K,K) = R^^{K) 

will be called an optimal perturbation of the method K for the linear problem. 

Taking K = 0 gives a (not perturbed) Runge-Kutta method ( |1.7| ) and a (not perturbed) 
stability function (j)K- In this case we denote the threshold factor i?Lin(Al, 0) simply by R{(j)K)- 
Clearly 

R{cfK) < R^^{K). (2.4) 

In the next section, we give upper bounds on R^^^(J^)- 

2.1. Upper bounds on the J;hreshold factor for optimal perturbations^. For any 

explicit perturbed RK method (AT, K) of linear order p, it can be seen that where 

fig p, with p < s, denotes the set of bivariate polynomials with the following properties: 

l.if(z,z)=J^^+ ^ cTjZ-^ + (z + z)'^(z,S); 
j=0 j=p+l 

2. dt is a polynomial of combined degree at most s — 1. 

In this section we investigate 

Rs,p = sup |i?(V') I ifiz, z) e ns,p| , 

Clearly 

R°^^{K)<Rs,p. (2.5) 


However, not all functions in H^^p can be realized as the stability function of an s-stage perturbed 
Runge-Kutta method (1.8), so inequality (2.5) is often strict (see Example 2.3 below). In case the 
optimal polynomial is realizable, the corresponding method may be of interest for the integration 
of linear systems. 


2 . 1.1 


we give an upper bound for Rg^p. In Subsection 


In Subsection ! 

to compute, Rg^p for given s and p, along with numerical values. 


2 . 1.2 


we give an algorithm 


2.1.1. Upper bound on Rg 


Lemma 2.1. Let (p(z) be a polynomial satisfying 

ip{z) = 1 7 iz H-h + 7 p+iz^’+^ H-h Jgz'" (2.6) 

Ij > 

j! 

Then the radius of absolute monotonicity of (p satisfies 

R{<p) < f/s{s-l)...{s-p+l). (2.7) 
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Proof. If R{ip) = 0, inequality (2.7) is trivial. Let ip{z) satisfy (2.6) and be absolutely monotonic 
at —r with r > 0. Then it can be written as 


^{z)= (i +fy=j2 £ 


1 t, 

Z J 


j=0 


j=0 \£=0 


\£ 


= E E- 

e=o \j=i 


where aj > 0, and = 1- As (/? is of the form (2.6), the coefficient of z'^ is larger than 1/p!. 

Some computations give 



Consequently, 


< ys(s- 1) • • • (s-p+1). 


We remark that equality in (2.7) is obtained for the polynomial 

where r = f/s{s — l)---(s—p+1). 

With the previous results, we can prove an upper bound on Rg^p. 
Theorem 2.2. 

Rs,p < ^s(s - 1) • • • (s - p + 1). 


( 2 . 8 ) 


(2.9) 


Proof If R{ip) = 0 for all then Rg^p = 0 and inequality (2.9) is true. Otherwise, there 

exists a function 'if G Ilg^p a.m. at (—r, —r) with r > 0. By [101 Lemmas 2.9 and 2.10], if is a.m. at 
the points (C,^), with ^ G [—^)0]- Writing 'ip{z,z) = and differentiating shows that 

all coefficients fijk are non-negative since ip is a.m. at (0,0). Thus ip{z,z) (viewed as a function of 
one variable) is of the form (2.6) and is a.m. at —r. Application of Lemma [2.1| gives the desired 
result. □ 

2.1.2. Numerical computation of bounds Rg^p. In this section we provide a means to 
compute tighter bounds on Rg^p for given s,p, using linear programming. The material in this 
section closely follows [151 Section 4.6.2]. The stability function (2.1) of an explicit s-stage perturbed 
Runge-Kutta method (1.8) with linear order p can also be written in the form 


S J 


j-0 i=0 




i + - 

r 


with jji = 


d^tpi 


j\ dz^~^dz^' 


( 2 . 10 ) 


furthermore, ip{z,—z) = 4>k{z) = exp(z) -I- 0{zP'^^). After considerable manipulation we find that 
ip{z, -z) = Yh=o where 


S j 

j—i m—ma,x{0,i—£) 


j-A/ £ \ (-l)- 

m ) \i — m) P 
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_ Table 2.1 

Rs,p: upper bounds on the threshold factors for optimal perturbations 



Hence we have the following problem for existence of a polynomial (2.10) with perturbed thresh¬ 
old factor at least r and order at least p: 


Given r find 7 such that 

> 0 0<i<j<s ( 2 . 11 a) 

Ci{r,^) = ^ 0<i<p. ( 2 . 11 b) 

i\ 


Since (2.11b) is a system of linear equations (in 7 ) then for any given value of r (2.11 1 represents a 
linear programming feasibility problem. Hence we can use bisection and an LP solver to find Rs,p, 
as was done for similar problems in [T51 [Ti] . Table 
up to ten. 


2.1 


gives computed values of Rg^p for s and p 


2.2. Examples. We now give some examples of optimal polynomials and Runge-Kutta meth¬ 
ods. 

2.2.1. Polynomials achieving i2s,p- The algorithm just described also provides coefficients 
for an optimal polynomial, which may or may not be realizable as the stability function of a 
perturbed Runge-Kutta method. 

The optimal Hrst-order polynomial for any s is just the stability polynomial of a (not perturbed) 
Runge-Kutta method consisting of s repeated forward Euler steps. 

The optimal order-two polynomial of degree s also has a simple form: 


^8.2(2:, 5 ) 


2{s + r) - 1 G ^ 2 
2(s -f r) V r 


1 

2(s -I- r) 



( 2 . 12 ) 


where r = Rs ^2 = Observe that bound (2.9) is sharp for p = 2. 

Some of the other optimal polynomials also have rational coefficients. Two optimal degree-four 
fourth order polynomials we found are 


tpl4iz,z) 
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and 


^Uz,'z)-l-^{l + l) +^(l + ^) (l + ^) +^(1 + 7) (l + ;)+^^ + 


where r = i? 4,4 = 2. Thus the optimal polynomial in Hg p is in general not unique. 

We have also computed (not shown in this paper) the exact values of Rs^p and polynomials i/'s.p 
for p = 3 and s = 3,4,5, 6 . In each case, we found an optimal polynomial of the form (2.10) where 
p of the coefficients jsj are non-zero and the coefficients ■jij with i < s are all zero. 

Exact values of Rs,p can be found in a systematic way as follows. First, a high-precision 
approximation can be computed using bisection and linear programming as described above. In 
practice, this yields a set of coefficients 7 ^; in which only p values are non-zero [7]. Setting the 
remaining values to zero a priori in (2.11b) yields a system of p -I- 1 equations in p unknowns which, 
nevertheless, possesses a solution. The solution may be found using a symbolic linear algebra 
package. We do not pursue this further in the present work. 



only upper bounds on what can he achieved. We do not pursue the topic further here. 


□ 


2.2.2. Optimal threshold factors for perturbations of specified 2-stage 2nd-order 
methods. We have no general method for finding nor a corresponding method. In this 

section we report results of some symbolic searches. In the case of the second-order methods, due 
to the small number of free parameters, it is not difficult to prove that the results below are truly 
optimal. 

Example 2.3. We consider explicit perturbed second-order 2-stage methods 


0 

0 

0 

a 

a 

0 

K 

1 - 2 S 

1 

2 a 


For these methods, function (2.1) can be expanded as 
, . . 1 , . . 


{K,K)\ 


where 


fill = b*Ae -\- FAe = 62021 + 62621 , 



0 

0 



621 

0 

(2.13) 

K 

61 

^2 

- z ) -\- fii{z -\- 

z) + / 32(2 + z)^ , 

(2.14) 

— bi 

^2 , 

(32 = b*Ae = 62621 

(2.15) 


The polynomial (2.14) is realizable (in the sense that it corresponds to a 2-stage Runge-Kutta 
method ( 2.13 [ ) ) if the first and last equations in (2.15) can he solved for 0,21 and 62 in K. A simple 
computation gives that a necessary condition is fifi — 2^2 > 0. Note that the stability function is 
independent of a. 

















With the help of Maple, we have computed the largest r such that (2.14) is a.m. at (— r, — r) 
and the polynomial is realizable. We have obtained that the optimal perturbation, denoted by K^, 
satisfies &2 = 021 = 0 and 61 = | (1/7 — 2 ), and 


^ (1 + ~ 1-21525 . (2.16) 

Observe that < i? 2,2 = '\/2. The stability function (2.1) for the optimal perturbed method 


- 9 (4 + v^) (1 + ^) +9(5-^) (1 + 


where r = 

Example 2.4. We consider perturbations of the classical four stage, order four method, of the 
form 


(2.17) 


0 

0 

0 

0 

0 


0 

0 

0 

0 

1 

2 

1 

2 

0 

0 

0 


0 

0 

0 

0 

1 

2 

0 

1 

2 

0 

0 


031 

0 

0 

0 

1 

0 

0 

1 

0 


041 

042 

0 

0 

K 

1 

6 

1 

3 

1 

3 

1 

6 

K 

bi 

^2 

0 

0 


We consider these perturbations because, in order to obtain a nonzero SSP coefficient for nonlinear 
problems, the analysis done in m shows that only the entries 031, 041, 042, 61 and 62 ^ K need 
be nonzero. To study SSP coefficients for the linear case, we have to analyze the perturbed stability 
function, that in this case is of the form 


\K,K) 


(z, z) — 1 + z + + -z^ + + z) + /?iiz(z + z) + /321 z‘^{z + z) 


(2.18) 


where 


/?! — + 62 j fill — g (^3 62 + 2 031 + 041 + 042 


). 


fi21 “ Y 2 ’ 


After some computations with (2.18), we obtain a coefficient « 1.66728 that is the 

positive root of the polynomial 15 a;"' — 4— 12 — 24a; — 24 = 0. The coefficients are 


fii = 


7r^ — 2r^ — 6r — 12 


fiii = 


5 — 2 r — 6 


fi21 = 


r — 1 


12 ’ " 12 ’ ' 6 
where r = RpifiK, K). With these values, the perturbed stability function can be written as 

701 (^1 + ^)+711(1+^) (^1 + ^) +721(1+^) (1 + ^)+740(1 + ^) , 

where 


7 oi = 


r{2 r^-r^- 6 ) 


6 


711 = 


(r2 + 2r-6) 


12 

9 


721 = 


r3 (j. _ 1 ) 

~6 


, 740 = 


24 















This perturbed stability function can be realized with the family of perturbations 


“31 = ^ (2 7’ - 2 - 042) , 

041 = ^ (5 - 6 r - 2 - 6 62) , 


3. Perturbed Runge Kutta methods for nonlinear problems. In this section we seek 
to answer the question posed in the introduction: for a given method K, what perturbation K gives 
the largest value of R{K, K)1 We begin by providing some upper bounds. 


It is convenient to write scheme (1.7) in canonical Shu-Osher form 


Y = VrUn + ar \ Y -\ - F 

r 


(3.1) 


where 


Vr = {I + rK) ^e, ar = r{I + rK) . 


(3.2) 


Observe that matrices K and ar have the same structure (strictly lower triangular, lower triangular 
or full). 

A Runge-Kutta method ([13 is said to be absolutely monotonic at r if (/ + rK) ^ exists and 
all entries of ar,Vr are non-negative. For a given method AT, the largest r such that the method 
is absolutely monotonic at r is known as the SSP coefficient, Kraaijevanger coefficient, or radius of 
absolute monotonicity m, and it is denoted herein by R{K): 

R{K) = sup {r I r = 0 or r > 0, (/ -I- rK)~^ exists, and > O} . (3.3) 


As usual, the inequalities above should be understood component-wise. 

Next, we consider perturbed Runge-Kutta methods (1.8). To study absolute monotonicity of 
perturbed Runge-Kutta methods, we write method (1.8) also in a canonical Shu-Osher-like form 


where 


Y — 


Y+-F 
r 



(3.4) 


'Jr = {I + rK + 2rK) ^e, (3.5a) 

= r{I + rK + 2rK)-^{K + K), (3.5b) 

Ao'^’^ = r{I + rK + 2rK)-^K. (3.5c) 


Observe that method (3.4), with 7 ^ = (I — aJIP — Q!(}°'^")e, is a perturbed Runge-Kutta scheme with 
Butcher coefficients 


K = -(I — al!P — a; 
r 


^)-i(a;iP-a^°™), K=-(/-a)lP-a; 

r 


down \ — 1 ^ down 


(3.6) 


10 












provided that (/ — q;“p — 0 ^!°™)“^ exists. 

The radius of absolute monotonicity of a perturbed Runge-Kutta method (RT, K) is the largest 
r such that 7 ^., and aj}®™ in (3.5) exist and are non-negative il Definition 3.1]: 

R{K, K) = sup |r I r = 0 or r > 0, (/ -I- rK 2rK)~^ exists, and 7 r, 0 ]!^, > o| . (3.7) 

3.1. Zero-well-defined perturbations. Regularity of (/ — aJfP — 0 :)!°™) is evidently impor¬ 
tant in our study. Observe that from (3.5) we have 


(/ - a"P - + rK) = (/ - 2a“™). 


(3.8) 


Consequently, if / -I- rK is regular for some r, then (I — a“P — af!°'"") is regular if and only if 
(/ - 2 a^°™) is regular. 

If / — 2a5(°'"" is singular, then the stage equations do not have a unique solution even for the 
trivial ODE given by / = 0. Hence we say that methods for which / — is singular are 

not zero-well-defined. See [5J Chap. 3] for the analogous definition in the context of traditional 
Runge-Kutta methods. 

3.2. Optimal perturbations. The optimal perturbed SSP coefficient of a Runge-Kutta method 
K is denoted by 


R°P^{K) = supR{K,K). 

K 


For a given method K that is (explicit/diagonally implicit/fully implicit), we consider the supremum 
over perturbations K that are zero-well-defined and correspond to the same class of methods. A 
matrix K such that R{K, K) = R°p^{K) is called an optimal perturbation. For some methods, the 
optimal perturbation is not unique. 

The following result shows that every method can be perturbed so as to give a method with 
strictly positive SSP coefficient. 

Theorem 3.1. Let K be a Runge-Kutta method that belongs to a specified class of methods 
(explicit, diagonally implicit, or fully implicit). Then it is always possible to find a perturbation K 
within the same class such that R{K,K) > 0. 

Proof. From [SJ Proposition 3.7], we have R{K,K) > 0 if and only if the Butcher coefficients 
satisfy 


K + K> 0, K>0, 


(3.9) 


and the following inequalities hold, 

Inc {{K + 2K){K + K)) < Inc {K + K), (3.10a) 

Inc {{K + 2K)k) < Inc {K ), (3.10b) 


where Inc (F) denotes the incidence matrix of matrix F defined as Inc (F) = (%, ) where qa = 1 if 

Consider first the implicit case. By making all entries of K positive we can satisfy (3.10), and 
by making them large enough we can satisfy (3.9). For the explicit and diagonally implicit cases. 
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note that if K,K are (strictly) lower-triangular, then the left-hand sides of (3.10) are also. Thus by 
making all the (strictly) lower-triangular entries of K positive, and by taking them large enough, 
we can satisfy the above inequalities. □ 

Observe that a perturbed Runge-Kutta method {K, K) can be interpreted as an additive 
Runge-Kutta method {K -|- K, K) for functions (/,/), and conditions ( |3.5| ) are the ones required 
for the absolute monotonicity of this additive scheme at (z, z) = (—r, —r) (see nni). From Lemma 
2.8 in [TU], we obtain that the stability function 4>(kk) defined by ( |2.1[ ), is absolutely monotonic 

(CiO = {—TT—f)- Consequently, 


R{K,K)<Ri^^{K,K)<R^Uk). 


(3.11) 


Furthermore, from SSP theory and inequality (2.4), we have 


R{K) < R{4>k) < R^liK), R{K) < R°^\K) < RI^Uk) . 


(3.12) 


The following example illustrates that R{(j>K) can be either larger or smaller than R°^^{K). 

Example 3.2. We consider the family of second order 2-stage Runge-Kutta methods (2.13) 
for a S K. For this family we have 


Wr > 0 
Or. > 0 


Thus 


a > 0 and 0 < r < — , 

a 


a > - and 0 < r < 


2a-1 


R(K) = < 


0 , 

2a- 1 
a 

1 

5 

a 


*/ a<2, 

if - < a < 1, 
if 1 < q; . 


(3.13) 


In Figure \3fl\ we show the threshold factor R((j)K) (thin solid blue line) and the SSP coefficient 
R{K) (thick solid black line). We also show the corresponding optimal coefficients for perturbed 
methods, namely, the optimal threshold factor R°^^J^K) (thin dashed blue line) and the optimal SSP 
eoefficient R°'^'^{K) for the perturbed method (thick dashed black line). 

We see that for optimal SSP method (a = 1) it is not possible to increase the SSP coefficient 
by means of perturbations. However, for a = (pjl — l)/2 it is possible to obtain a perturbation that 
raises the SSP coefficient to R°^^(K) = R°iiJyK) = (l -|- -s/?) ~ 1.21525 (see (2.16)/ 

We also see that for 2/3 < a < 1 we obtain R{(pj^) < R°'p^{K), whereas for 0 < a < 2/3 and 
for I < a we have R°^^{K) < R(4 >k). 

Coefficients of the perturbations that give rise to these values are given in Appendix\5 . 1\ 


3.3. Upper bounds on the SSP coefficient for perturbed RK methods. In this section, 
we explore some upper bounds on the SSP coefficient R°p^{K) where K is an s-stage order p method. 
A straightforward upper bound is obtained from inequality (|3.12 ) and Theorem 2.9 


r°pCK] < ifs 
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Fig. 3.1. Family of second order 2-stage methods: SSP coefficients for unperturbed methods and optimal SSP 
coefficients for perturbed methods. 


As the next Theorem shows, the largest positive value such that vector Vr in (3.21 is non-negative 
is also an upper bound for R°p*'{K). 

Theorem 3.3. Consider an explicit Runge-Kutta method K and let Vf, be the largest positive 
value such that vector Vr in (3.2) is non-negative. Then 


R°'P\K) < re. 


(3.15) 


Proof. Let r = R°^^{K). Then 7^ = (/ — > 0, and thus from (3.2) and (3.8) we 


get 


(/ - > 0 . 


(3.16) 


As > 0, and since we consider only explicit, zero-well-defined perturbations, I — 20)?°'^" is 

an M matrix. Thus (/ — 2 q!(}°™)“^ > 0. If we multiply (3.16) by (/ — 2 q;(}°™)“^ we obtain that 

Vr >0. □ 

From Theorem [3^ we obtain that 


R{K) < R°P\K) < re 


(3.17) 


Consequently, for those methods such that R{K) = rg, the SSP coefficient cannot be increased 
by perturbation. This is the case for the family of second-order two-stage methods. For a > 1, 
R{K) = re = 1/a (see Example |3.2| ). 

On the other hand, if R{K) < Ve one can try to find a perturbation to increase the SSP 
coefficient. This is the case for the classical 4-stage order 4 method for which R{K) = 0 and 
re ~ 1.2956, the real root of — 2x^ -I- 4a; — 4 = 0. 

Another interesting bound, for explicit methods only, can be obtained in terms of the Butcher 
coefficients of the Runge-Kutta method K. 
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Theorem 3.4. Consider an explicit Runge-Kutta method {K, K) with perturbed SSP coefficient 
R°P^{K) >0. Let K = (a^). Then 


R°P\K) < 


maxii \ai 


(3.18) 


Proof. The proof is similar to that of Lemma 3.2]. Consider an optimal perturbation K and 
set r = R°P^{K, K) > 0; consider too the canonical representation (3.4). Let A = aJlP + = 


(a.,), r = aJlP/r = L = 


,,down 


/r = (fdij)] observe that A,r,r > 0, and that A = r(r + f). 


As (/ — A)e = 7r > 0 and atk > 0, we have aik < 1; as (/ — A)K = L — L, we have 


— (dik ffik T ^ ^ OLijOjk 

j=k+l 


(3.19) 


As aik = r{j3ik + /3zfe), then jd^k + fdik = aik/r < l/r. In particular, from (3.19), 


|a2i| = 


/321 — ,021 


< (d2i + id2i < - ■ 


We proceed by induction on row i of K. Assume that ja^ | < l/r, for * = 2, ...,£, j = 1, ...,£— 1, 
and consider row I +1. Then, from (3.19), 


Wi+1,11 


I 

/3^+i,i ~ Q!^+i, i Oj,! 

1=2 


I 

< ldf.+i,i + Pt.+i,i + o:i+i,j \ Oj,! I 
1=2 


< 


-a£+l 


1 ^ 


a^+1,1 


1=2 



£ 

E “^+1 

1=1 



A similar argument can be used to show that |a^+ij | < l/r, j = 2,..., ^. The Theorem follows by 
induction. □ 

Consequently, 


R{K) < R°^\K) < -^p. 

maxy \aij\ 


For those methods such that R{K) = l/max^ \aij\ it is not possible to increase the SSP coefficient 
by perturbing the method. This is the case for all known optimal explicit SSP RK methods of 
orders one through four, with any number of stages m- 

For the restricted class of perturbations considered in [22] , similar results were obtained in [22] 
Theorems 3.1, 3.4, 3.5 and 3.6]. Theorem 3.4 extends those results, showing that no improvement 
in the radius of absolute monotonicity is possible for many optimal SSP methods, even when more 
general perturbations are considered. 


3.4. Relations among the Butcher and canonical Shu-Osher representations. For 

perturbations of explicit Runge-Kutta methods (or of other methods with one stage equal to Un) 
there exists a certain simple transformation that may yield a larger value of R{K, K) - and that 
never yields a smaller value. 
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Proposition 3.5. Let an s-stage explicit perturbed Runge-Kutta method be given with coeffi¬ 
cients 7r, aJJP, > 0, where r is the radius of absolute monotonicity of the method. Consider 

the perturbed method with coefficients 


7 = ( 1 , 0 ,..., 0 )*, 
s-? = <? + (>y2 

gdown ^ ^down ^ 


(3.20a) 

2<i<s (3.20b) 

2<i<s. (3.20c) 


Then the perturbed method with coefficients {'jr, a(!°™) and the modified perturbed method with 

coefficients (7, a“P,y°™) correspond to the same RK method K. The modified perturbed method 
has radius of absolute monotonicity at least equal to r. 

Proof. It is easily seen that the modihed method is equivalent to the original one when f = f, 
so they correspond to the same unperturbed method. Meanwhile, the transformation never leads 
to negative coefficients, so the modified method is a.m. at r. □ 

Remark 2. Proposition \3.^ is also valid for Runge-Kutta methods whose first row is equal to 
zero. ^ 

In the Butcher form (1.8) it is obvious which perturbed methods {K,K) correspond to a given 
method K. In the canonical Shu-Osher form it is less obvious. The following lemma characterizes 
which methods of the form (3.4) are perturbations of a given method (3.1). 

Lemma 3.6. If method (3.4) is a perturbation of method (3.1), then their coefficients are related 
as follows: 


(/ - 2ay")a, = (a"P - a^°™) 

(/-2y°™)n, = 7.. 


(3.21a) 

(3.21b) 


Furthermore, if (3.21) holds and the perturbation is zero-well-defined, then (3.4) is a perturbation 

of (ra. 


Proof To prove the first part, take / = / in (3.4) to obtain: 

y = 7,u„ + (a;iP + a^°™)y + (a“P - 0^")-^. 

r 

Subtract 2a)J°™y from both sides to get 

(/ - 2ad°™)y = jrUn + (ar - ar°™) • 


(3.22) 


Substituting (3.1) in the above gives 


{I - 2ay")z;,u„ + (/ - 2a“™)a, [Y +-F ] = 7,n„ + (a"P - a 


down 


)lr + y 


Equating coefficients yields (3.21). 


To prove the second part, assume I — 2a)!°™ is invertible and write (3.211 as 

ar = {I- 2y°™)-i«p - y°™) 
i;, = (/-2y°™)-V,. 


(3.23a) 

(3.23b) 
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Substitute (3.23) in (3.1), multiply on the left by (/ — and follow the steps above in 

reverse. □ 


Lemma 3.6 does not imply that the perturbation (3.4) is unique for a given r; see Proposition 


Remark 3. The necessity of the zero-well-defined condition in the second part of Proposition 
can be seen from the following example. We take the implicit trapezoidal Runge-Kutta method 


0 

0 0 

1/2 

1/2 1/2 


1/2 1/2 


The canonical form (3.1) is then 


j 

( 0 

0 


1 j 

( ^ \ 

Or = 

r 

r+2 

r 

r+2 

0 ' 

Vr=\ 

2-r 1 

1 

1 

r 

o) 

> 1 

1 2-r j 


\r+2 

r+2 


\r+2/ 


Then (3.21) is satisfied - for any r - by 



/1/3 0 0^ 

aup = ad°wn= 0 1/2 0 

\ 0 1/2 Oy 

However, this method - which involves a perturbation that is not zero-well-defined - is not a per¬ 
turbation of the original method. □ 

3.5. Computing optimal perturbations. In this section we present two algorithms to sym¬ 
bolically or numerically find the optimal perturbed SSP coefficient and a corresponding perturbation 
of a given RK method. The first algorithm is proven to approximate the optimal value to any accu¬ 
racy, contingent on the computational solution of linear program subproblems. It is only valid for 
explicit perturbations. The second algorithm is analytical and exact, and valid for both explicit and 
implicit methods, but it is not proven to give the optimal value. The results of the two algorithms 
coincide (to high precision) for all explicit methods on which we have tested them. 

3.5.1. Provably optimal algorithm for explicit perturbations. In the foregoing, we 
have shown that Hnding an optimal perturbation consists of determining the largest r such that 


there exists a splitting satisfying (3.21) with positive coefhcients. Note that the range of values 


for which a method {K,K) is absolutely monotonic is always the interval [0, R{K, K)]. Therefore, 
one way to find the largest r is to devise a method for testing for a given r whether there exists 


a perturbation K such that R(K,K) > r. For given method (1.7) and value of r, the system of 
equations (3.21) together with the inequalities > 0 constitutes a linear programming 


(LP) feasibility problem. The following theorem is an immediate consequence of Lemma 3.6 


Theorem 3.7. Let an s-stage RK method K and a positive number r be given. There exists a 
perturbation K with R{K, K) > r if and only if there exists an (s -I- 1) x (s -I- 1) matrix such 

that (/ — is regular and the following componentwise inequalities hold: 


I ^ down ^ n 

+ > L) 


(/ - 2af 

{I - 2af°'^^)vr > 0 


down \ n 
> 0 . 


(3.24a) 

(3.24b) 

(3.24c) 
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The linear program (3.24) can be solved by standard LP solvers. By embedding this solution in a 
one-dimensional root-finding algorithm, optimal perturbations can be found. An algorithm based 
on bisection follows. 


Algorithm 1 Optimal explicit perturbation 


Input: K 

^max ■— 1/ laiax II , Tniin ■— 0- 


while 


^ m a.v ^ n 


> e do 


rp _ max ~l~^min 

' ~ 2 


Compute the coefficient matrices ar,Vr using (3.2) 


Solve the LP given by (3.24). 
if it is feasible then 


else 




end if 
end while 
return 


Assuming the solution of the LP is correct, the algorithm provably finds an optimal explicit 
perturbation. However^for implicit perturbations the LP solver may converge to a solution (like 
the method in Remark ^ above) for which I — is singular. 

3.5.2. Iterated splitting algorithm. We next investigate how to choose directly 

so as to find a perturbation with radius of a.m. at least r. The following result suggests an approach. 


Lemma 3.8. Given an explicit Runge-Kutta method (3.1), let aJfP > 0,0))°™ > 0 denote 
coefficients of a zero-well-defined perturbation of (3.1). Then there exist matrices > 0, a~ > 0 
such that 


0 “'' = (/ -t- 2a ) *a 


^down 


= {I + 2a-) 




and ar = a~^ — a . 

Proof. Since the perturbation is zero-well-defined, we can define 


= 

oT = {I — 2a 


a ■ = (7 - 2a^°™) *a. 


i\ —i dc 
) “r 


down\—1 up 
) 

down \ — 1 _ down 


(3.25a) 

(3.25b) 


(3.26a) 

(3.26b) 


Then, by (3.21a), ar = a’*' — a . Furthermore, since I — 2a))°™ is an M-matrix, we have a’*' > 0 
and a~ > 0. Solving (3.26) for a“P,a((°'^" gives (3.25). □ 


In the next algorithm we use the following notation: 


0 


if Xij > 0 
if Xij < 0. 



if Xij > 0 
if Xij < 0 , 


and thus x = (a;)"*" — (x) is a sign splitting of matrix x, with (x)"*" > 0, (x) > 0. 
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Given a perturbed Runge-Kutta method (3.4| with 7^ = ei, where ei = (1,0,... ,0)*, and 


or containing negative values, we construct 

> = (/+ 2 «P)-+ 2 (a^own)-)-! 


7 = (/ + 2 (a^)- + 2 (a^°™)-) ^ («p)+ + (a^°™)-) 


a 

-down 


-down = (J + 2 (a^)- + 2 (a^°™)-) ^ ((a"P)- + (a^°™)+) 


where a“P = (aJlP)^ — 
exists. Using Lemma 


K!P)- 


.,down 


3.6 




,,down\ — 


) , provided that / + 2(a^P) 


(3.27a) 

(3.27b) 

(3.27c) 

down\ — 


it is straightforward to prove that, if method q;“p, is a perturbation 

of (|3.l|), then (|3.27|) is also perturbation of (|3.1|). Next, for explicit methods, we perform trans¬ 


formation (3.20). In this way, (3.27) followed by transformation (3.20) give a perturbation of the 
form (3.4) with 7^ = 61, that we denote by d“P, If d“P > 0, > 0, then r is an SSP 


coefficient; otherwise, we can repeat the above process. 

The following lemma studies the sign of (dJlP)ij, (d(!°'^")y when (aJlP)^ < 0 or (0)}°™)^^ < 0. 
For the sake of clarity, we drop the index r. 

Lemma 3.9. We consider a perturbed explicit Runge-Kutta method with coefficients 7 = ei, 


q;“p, and the perturbation d“P, obtained by computing (3.27) followed by transforma¬ 


tion (3.20). Assume that jo > 2 is the first row with negative terms in q;“p or 
the largest index mo > 1 such that ct^^rno < ^ or < 0. Then 


Let mo be 


1. For first to {jo — l)-th row, we have: qAj = and = 0 for 1 < i < jo — 


* j 


1 < J < Jo - 2. 

2. For the jo-th row, we have: 

(a) If mo = 1, then < 0 or < 0. 

(b) If mo > 2, then, > 0 and <5;^°™ > 0. 

(c) For 1 < mo < jo — 2, we have > 0 and > 0 for £ = mo -\-1, ■ ■., jo — I- 

Proof. If jo is the first row with negative terms in a“P or straightforward computations 

give that d“j = a'fj and = 0 for 1 < i < jo — 1, 1 < j < jo — 2, and 


Jo-1 


“lb = “Zi - 2 E («.)- + «T")b <!. 


(3,28a) 


io-1 


^Zi = (^Z,iZ + «r)- - 2 E azzy + «rr) ^=2,...,jo -1. (3.28b) 

i=,^+l 


and 


down ^ down 
« J 0.1 =«^ 0.1 


Jo-1 

2E((“Z.)^ + «T)b 


i-2 


^ down 

o:^,l , 


down 


jo-1 

E 

i=e-hi 


(«Z<r + «r)*-2 E («,.)■ + «”")■) 


^ down 


(3.29a) 

£ = 2,...,jo-l. 

(3.29b) 


Let mo be the largest index mo > 1 such that < 0 or 0;^°™^ < 0. In this case, q;“U > 0, 
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> 0 for i = mg + 1,..., jo — 1, and thus 


“in,7 ’ ) 


down\+ _ down 


Qio 


K",7)- = «°rr = o> 


j = mo + 1 ,..., jo - 1 ■ 



If mo = 1, from p.28a| ) and p.29a[ ) we get = ctj^i and and thus &jgi < 0 or 

"down , 

“lO,l 


5,up _ 
jo,mo 


= r + I” > 0 

^ Jo,rngj ^y'^Jo,mo) — 


"down ^ / up 
“lo.mo '.'-^10 


= )- + )+ > 0 


Finally, for 1 < mo < jo — 2, from (3.28b) and (3.29b I we get that, for £ = mo + 1,..., jo — 1, we 
have 


- up 

^0 ^ " 


= «,er > 0 


down /^down\+ n 

^jo,e = i^jo,i ) ^ 0 


□ 

Consequently, if matrices a“P and contain negative elements in the second or later 

columns, an iterated construction of perturbations d“P, removes these negative values ob¬ 

taining a perturbation with non-negative elements from second column on. However, if in a row jo 
we have: 


<r.i<o 

and 

a;P,>0 £ = 2 ,...,jo-l. 

(3.30) 

or 




^ down ^ n 
«10.1 <0 

and 

<7>0 ^ = 2,...,jo-l. 

(3.31) 

the new perturbation dJlP, q:^°™ 

will also contain negative elements in the first column. 



We now give Algorithm to determine whether there exists a perturbation with a.m. radius r 
for a given method. 


Algorithm 2 Existence of a perturbation with radius r 


Input: r, K 

Compute the coefficient matrices ar,Vr using (3.2). 

Set q;“p = Or and 0 '^°™ = 0. 

while a'^P or has any negative entries do 


If K has a zero row, perform the transformation (3.20). 

If a“P, Q-down > stop. This is a feasible perturbation. 

If condition p.30|) or (3.31) hold, stop. A feasible perturbation cannot be found. 
Set a- = (a“P)- + (a^^+ and a+ = (a"P)+ + ( 0 '^°™)- 
Compute a new splitting: 


a"P = (/-t2((a"P)-^ a+ 
adown = (/ + 2((a"P)- -h a" 


end while 
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Order 

Stages 

Method 

R{K) 


Bound (3.18) 

Bound (3.14 

) Property C 

1 

1 

Forward Euler 

1 

1 

1 

1 

True 

2 

2 

Midpoint 

0 

0.732 

1 

1.414 

True 


2 

Minimal trunc. error 

0.5 

1 

1.333 

1.414 

True 


2 

SSP22 ng 

1 

1 

1 

1.414 

True 


2 

SSP22* [g 

0.784 

1.215 

1.215 

1.414 

True 

3 

3 

Heun33 [8] 

0 

0.776 

1.333 

1.817 

False 


3 

SSP33 |23] 

1 

1 

1 

1.817 

True 

4 

4 

RK44 (Kutta) 

0 

0.685 

1 

2.213 

False 


5 

Merson (TH] 

0 

0.242 

0.5 

3.309 

False 


10 

SSP104 [ig 

6 

6 

6 

8.425 

False 

5 

6 

Fehlberg [3] 

0 

0.057 

0.125 

3.727 

False 


7 

Dormand-Prince [g 

0 

0.040 

0.086 

4.789 

False 


8 

Bogacki [T] 

0 

0.313 

0.859 

5.827 

False 


7 

SSP75 ng 

0 

1.396 

1.792 

4.789 

False 


8 

SSP85 |2g 

0 

1.875 

1.919 

5.827 

True 


9 

SSP95 |2g 

0 

2.738 

3.198 

6.853 

False 

6 

9 

Calvo [g 

0 

0.021 

0.059 

6.265 

False 

8 

13 

Prince-Dormand [2g 

0 

0.013 

0.059 

9.212 

False 


Table 3.1 

Properties of some RK methods and their optimal perturbations. The optimal perturbed radius of absolute 
monotonicity was computed by both the linear programming algorithm and the iterated splitting algorithm; in every 
case they gave identical results (up to roundoff errors). Decimal values have been truncated to the number of digits 
shown. 


Remark 4. This approach seems to lead to optimal splittings for all the explicit methods on 
which we have tested it. However, for all implicit methods we have tested, it fails to increase the 
radius of absolute monotonicity at all. Even for explicit methods, we have no proof that it’s optimal 
because one could use + S, {ar)~ + 6, in place of {ur)'^, where 6 is any non-negative 

matrix. □ 


3.6. Examples. In this section we compute optimal perturbations of some existing methods, 
using the algorithms described in the last section. 

We have computed optimal perturbations for several known explicit methods using the two 
algorithms described above. In all cases, the two algorithms gave the same values. It thus seems 
possible that Algorithm 2 also gives truly optimal results in general, but we do not have a proof. 
Properties of the methods studied are given in Table [O] Several interesting facts are evident: 

• For all optimal SSP methods (up to order four), perturbation cannot yield a larger co¬ 
efficient. This is evident already from the bound (3.18). For all other methods, some 
improvement is achieved. 

• Consistent with Theorem 3.4 for every method considered, it is possible to achieve i?°P* > 0 
by some perturbation. 

• The simple bound (3.18) predicts the optimal coefficient to within a factor of three in every 
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case. 

• The methods SSP75, SSP85, and SSP95 are optimal methods found in [52], with property 
C. By considering methods without property C, we obtain slightly larger coefficients for 
perturbations of SSP75 and SSP95. On the other hand, relaxing the column assumption 
gives no benefit in the case of the SSP85 method. 

• The values found have been truncated to three decimal places but are known to greater 
precistion. For the 4-stage, order-four method of Kutta, the three-digit value of R°p^{K) 
given in the table matches the value found by Shu and Osher. However, the exact (irra¬ 
tional) value is slightly larger and is given in the appendix. 

4. Conclusions. In this work we have studied SSP coefficients for perturbations of a given 
explicit Runge-Kutta method. We have considered both the linear and the nonlinear case, and 
have obtained useful bounds on the threshold factor and on the radius of absolute monotonicity 
for perturbed Runge-Kutta methods. We have also provided an algorithm for computing optimal 
perturbations of explicit Runge-Kutta methods, and given optimal perturbations for many methods 
from the literature. 

This work seems to provide a complete picture for the case of most interest: explicit methods 
applied to nonlinear problems. Nevertheless, some other interesting issues remain unsolved. These 
include: 

• A method to compute optimal perturbations for linear problems. 

• An algorithm for obtaining optimal splittings of implicit methods. 

These may be a starting point for future work. 


5. Appendix. In this section we give additional details on SSP coefficients and optimal per¬ 
turbations of second order 2-stage Runge-Kutta methods and the classical 4-stage fourth order 
Runge-Kutta method. 


5.1. Second order 2-stage methods. We consider the family of 2-stage second order meth¬ 
ods (2.131. In example 2.3 we studied perturbations that increase the SSP coefficient for the linear 
case. For nonlinear problems, in example |3.2[ figure [ tT] shows the values of for a € [—3, 3]. 

In this section, for each a, we give the expressions for R°p*^{K) and we show optimal pertur¬ 
bations KjsfL such that R{K,Kjqi,) = R°^^{K). It is important to point out the convenience of 
choosing = Kl, where Kl denotes the optimal perturbation for the linear case. In this case, 
we have not only R{K,Kl) = but also Ri^[a{K, Knl) = The computations re¬ 

quired to obtain the results in this section have been done with the symbolic computation program 
Mathematica. 

If we denote by r = i?°P*(Ar), we have that 



if 





(5.1) 

Next we give optimal perturbations A^vl- 

For a < 0, we obtain that it is not possible to obtain a perturbation of the form ( |2.13| ) with 
62 = 0 and 021 = 0. Consequently, ^ and we always have that i?Liii(A, Ajvl) < 


— 1 -(- Oi -(- '\/3cr^'^^'2c(~-t-^ 


if a S 


a 
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Optimal perturbations of the form (2.13) for different values of a < 0 must satisfy the following 
conditions. ^ 

• For —5 (l + V7) < a < 0, the coefficients 021 , bi and 62 in Knl must satisfy 


—a < 0,21 < 


1 — r a 


2 r 


&i =- 


ra2i 

2a 


62 — — 


2 a ’ 


where r = R°p\K). 

For a < — 5 (1 + a/T), we should have 


1 ~ -2a2-2a + l 

021 =-a, -7^<bi< --- 

2 a 4 a 


1 ~ 2abi-l 

-< 62 < - - - 

2 a 4 a 


For a > 0 we can find optimal perturbations with &2 = 0 and 021 = 0. Coefhcient bi must satisfy 
the following conditions. 

• For 0 < a < (-1 + Vf) /2, we have that 


bi = 


V3a2 - 2 a + 1 - 


a 


2a 


(5.2) 


Thus there is a unique of the form (2.13). In this case, we have R{K) < R{K, = 

• For (-1 + V7) /2 < a < 1, we also get R(K) < but in this case the optimal 

perturbation Knl fs not unique. All the perturbations with bi satisfying 

1 — a ~ 2 a^ — 2 a + 1 

- < 61 < , 

a 4a 


are optimal. In particular, we can take ATatl = Kl- With this choice, R{K,Kl) = 
R°P^K) = Ija and i?Lin(A:,^L) = « 1.22. Furthermore, a = (-1 + 7?) /2 

provides the largest SSP coefficient within the family of 2-stage second order method (see 
figure 3.1). 

For 1 < a, we have R{K) = R°p*^(K) = 1/a and the optimal perturbation ATatl is not 
unique. All the values 


0 < 61 < 


2a^ -2a + l 
4a 


give optimal perturbations. We can take ATatl =0, but in this case i?Lin(Af, 0) < (AT). 

A better choice is K^l = K^. Observe that, for a = 1, we get the optimal SSP coefficient 
R{K) = 1 that cannot be increased by perturbations. 

Next, we consider some concrete values of a to show the the expressions of the perturbations. 
For each value, we give the Butcher tableau of the perturbation and matrices a“P and in 

([ 33 . 

• For a = 1/2 we get method RK2a in [12] with R{K) = 0. With perturbation 


^/00 0\ /00 0\ /O 0 0\ 

iF = 0 0 0 , a"P = 5i 0 0 , a‘^°™ = 0 0 0 ,7= l-&i , 

\bi 0 0/ Vo 2 bi 0/ Vl-261 0 0/ \ 0 / 

where 61 = 5 (a/S — l), we get R{K, K) — i?Ljn(A', K) = a/S — 1. 


22 















• For a = 2/3, we have a nontrivial SSP coefficient = 1/2, but we can increase this 

value to R{K, Ki) = R°p^{K) = 1 with perturbation 



0 

0 \ 

/ 0 

0 

o\ 


Ki=\ Q 

0 

0 

, a“P = 2/3 

0 

0 

down 

, a 

Vl/4 

0 

0 / 

V 0 

3/4 

o; 



/ 0 0 0 \ 
0 0 0 
\l/4 0 0/ 



For this perturbation, R{(f)x) 
the first column of and a 


according to (3.20), 


7 


( 1 , 0 , 0 )* by modifying 


_ / 0 0 0 \ / 0 0 0 \ / 0 0 0 \ /l\ 

if 2 = 1/6 0 0 ,a"P= 5/6 0 0 , = 1/6 0 0 , 7 = 0 

\3/8 0 0/ \ 0 3/4 0/ \l/4 0 0/ \0/ 


• As it has been pointed out above, the largest value in the a-family of 2-stage second order 
schemes is R°p^{K) = (1 -b ■\/7)/3 and it is obtained for a = (V7 — l)/2. The perturbation 
is of the form (2.13) with 61 = (V? — 2) /2, and 


/O 0 0\ 

1 0 0 , 

Vo 1{4 + V7) Oj 


/ 0 0 0 \ 
0 0 0 
Vi(5-V7) 0 oJ 



This is the perturbation obtained in [6l Table V] by numerical search in the class of per¬ 
turbations considered in [5]. 

5.2. Classical fourth order 4-stage method. For nonlinear problems, applying the analysis 
above, we find that the optimal perturbation of the classical method has SSP coefficient given by the 
real root of + 2x^ -I- 4a: — 4 = 0, which is approximately R°p^{K) « 0.685016. The corresponding 
perturbation is not unique. For instance, we can take 7 ^. = (1,0, 0,0,0), and all entries of 
equal to zero except 


„2 2 

/ down\ ' / down\ ' /r o\ 

[ar )31 = ^ (ttr )42 = y> (5-3) 

where r = i?°P‘(Ar). However, there exist other optimal perturbations with additionally (q :/°™)42 = 
e where 0 < e < 0.782. 

We remark that nearly-optimal perturbations for this method are given in [231 P- 448] and [llj . 
Interestingly, these different perturbed methods have different values of i?Lin(7fj Tf). 
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